#include "AnaDef.h"

const Double_t EPSILON = 1e-10;

/*
//cutElectron_MCjet_19012014
Hh THnIDFF content:
THnSparseT<TArrayD> (*0x2c621b0): "THnIDFF" ""
  11 dimensions, 8.9389e+06 entries in 3225661 filled bins
THnIDFF axis 0: TrackEta, nbin= 4
THnIDFF axis 1: JetPt, nbin= 20
THnIDFF axis 2: TrackTPCsig, nbin= 1200
THnIDFF axis 3: Log10TrackP, nbin= 300
THnIDFF axis 4: Log10TrackPt, nbin= 50
THnIDFF axis 5: z, nbin= 25
THnIDFF axis 6: xi, nbin= 50
THnIDFF axis 7: pdg, nbin= 7
THnIDFF axis 8: comb, nbin= 7
THnIDFF axis 9: tof, nbin= 7
THnIDFF axis 10: tpc, nbin= 7

 */

const Int_t axEta = 0;
const Int_t axJpt = 1;
const Int_t axTPC = 2;
const Int_t axMom = 3;
const Int_t axPt  = 4;
const Int_t axZ   = 5;
const Int_t axXi  = 6;
const Int_t axPDG  = 7;
const Int_t axComb = 8;
const Int_t axtpcid = 10;
//actually this is "comb"
Int_t axPreID = -999;

const Int_t MCoffset = 4;

TH1D * GetHyield(TH2D *h0)
{
  //only consider the dEdx in range
  return h0->ProjectionX(Form("%syield", h0->GetName()), 1, h0->GetNbinsY(),"oe");
}

void GetHist(const TList *ll, THnSparseD * &hh, const TString name)
{
  hh=(THnSparseD*)ll->FindObject(name);
  if(!hh){
    printf("no hh %s!!\n",name.Data()); exit(1);
  }
  printf("\nHh %s content:\n", name.Data());
  hh->Print();
  const Int_t ndim = hh->GetNdimensions();
  for(Int_t ii=0; ii<ndim; ii++){
    const TAxis *ax = hh->GetAxis(ii);
    const Int_t nbins = ax->GetNbins();
    printf("%s axis %d: %s, nbin= %d\n", name.Data(), ii, ax->GetTitle(), nbins);
  }
  printf("\n");
}

void ResetMC(THnSparseD * hh)
{
  hh->GetAxis(axPDG)->SetRange(MCoffset,1000);
  printf("ResetMC setting MC range %d %d, %e %e\n", hh->GetAxis(axPDG)->GetFirst(), hh->GetAxis(axPDG)->GetLast(), hh->GetAxis(axPDG)->GetBinLowEdge(hh->GetAxis(axPDG)->GetFirst()), hh->GetAxis(axPDG)->GetBinUpEdge(hh->GetAxis(axPDG)->GetLast()));
}

void OutputHist(THnSparseD * hh, const TString tag)
{
  TList *ll=new TList;
  ll->SetOwner(kTRUE);

  //===========================================

  TH2D *kk = hh->Projection(axPt, axMom); ll->Add(kk);
  kk->SetName("kk");

  TH2D * zz = hh->Projection(axZ, axMom); ll->Add(zz);
  zz->SetName("zz");

  TH2D * xi = hh->Projection(axXi, axMom); ll->Add(xi);
  xi->SetName("xi");

  if(tag.Contains("kMC1")){
    TH2D * mcpt = hh->Projection(axPDG, axPt); ll->Add(mcpt);
    mcpt->SetName("mcpt");

    TH2D * mczz = hh->Projection(axPDG, axZ); ll->Add(mczz);
    mczz->SetName("mcz");

    TH2D * mcxi = hh->Projection(axPDG, axXi); ll->Add(mcxi);
    mcxi->SetName("mcxi"); 
  }
  //===========================================
  const Int_t nt  = axPreID>=0 ? hh->GetAxis(axPreID)->GetNbins() : -999;
  const Int_t tup = axPreID>=0 ? hh->GetAxis(axPreID)->FindBin(-1): -999;
  printf("OutputHist hh %s tag %s nt %d tup %d\n", hh->GetName(), tag.Data(), nt, tup);

  if(nt>0){
    //TOF -1
    hh->GetAxis(axPreID)->SetRange(1, tup);
  }

  TH2D * horiginal = hh->Projection(axTPC, axMom); ll->Add(horiginal);
  horiginal->SetName("horiginal");

  TH1D *horiginalyields = GetHyield(horiginal); ll->Add(horiginalyields);

  if(tag.Contains("kMC1")){
    TH2D * mcmom = hh->Projection(axPDG, axMom);
    ll->Add(mcmom);
    mcmom->SetName("mcmom");


    for(Int_t itype=0; itype<4; itype++){
      hh->GetAxis(axPDG)->SetRange(MCoffset+itype, MCoffset+itype);
      TH2D * tpcmc = hh->Projection(axTPC, axMom); ll->Add(tpcmc);
      tpcmc->SetName(Form("tpcmc%d", itype));
    }
    ResetMC(hh);
  }

  if(nt>0){
    //TOF 0 1 2 3
    hh->GetAxis(axPreID)->SetRange(tup+1, nt);

    TH2D * tof = hh->Projection(axPreID, axMom); 
    ll->Add(tof);
    tof->SetName("tof"); 

    if(tag.Contains("kMC1")){
      TH2D * mctof = hh->Projection(axPreID, axPDG);
      ll->Add(mctof);
      mctof->SetName("mctof");
    }

    //TOF any
    hh->GetAxis(axPreID)->SetRange(1, nt);
  }

  
  //===========================================

  const Double_t thres = 1000;
  if(horiginal->GetEntries()>thres){
    printf("OutputHist hh %s tag %s nn %.0f\n", hh->GetName(), tag.Data(),horiginal->GetEntries());

    const TString dirname(Form("histRec%s_LHC%s", hh->GetName(), tag.Data()));
    gSystem->Exec(Form("mkdir %s", dirname.Data()));
    TFile * fout = new TFile((dirname+"/")+(dirname+".root"),"recreate");

    ll->Write();
    
    fout->Save();
    fout->Close();
    delete fout;
  }

  //protected by TH1::AddDirectory(kFALSE);
  delete ll;
}

void getOutputHist(THnSparseD * hh, const TString tag)
{
  printf("tag: %s\n", tag.Data());

  //====================================================
  if(tag.Contains("kMC1")){
    ResetMC(hh);
  }
  //====================================================

  const Int_t nEta = hh->GetAxis(axEta)->GetNbins();
  const Int_t nJpt = hh->GetAxis(axJpt)->GetNbins(); 

  //should not be any over-/under-flow!!
  TH2D * hall = hh->Projection(axEta, axJpt);
  for(Int_t ieta=0; ieta<=nEta+1; ieta++){
    for(Int_t ijpt=0; ijpt<=nJpt+1; ijpt++){
      const Double_t icont = hall->GetBinContent(ijpt, ieta);
      printf("%s Entries eta %d, jpt %d: %.0f\n", hh->GetName(), ieta, ijpt, icont);
      if(ieta==0 || ieta==nEta+1 || ijpt==0 || ijpt==nJpt+1){
        if(icont>EPSILON){
          printf("over-underflow non-0!! %s %s %d %d %e\n", hh->GetName(), tag.Data(), ieta, ijpt, icont);
        }
      }
    }
  }

  TFile * fall = new TFile(Form("%s_LHC%s_allEta_allJpt.root", hh->GetName(), tag.Data()),"recreate");
  hall->Write();
  fall->Save();
  fall->Close();
  delete fall;
  delete hall;


  //===========
  TAxis *jax = hh->GetAxis(axJpt);

  const TString hn = hh->GetName();

  //const Int_t firstjpt = hn.Contains("IDFF") ? 2 : 0;

  //test fine binning in 5-20 in 15 bins
  //const Int_t jptbin0[]={1,6, 11,  1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
  //const Int_t jptbin1[]={5,10,15,  1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};

  //9-10 GeV/c ETA2 still enough stat. than 15-20 eta2
  const Int_t jptbin0[]={1,2,3,4,5, 1,1, 1,   6, 6, 6, 6, 11};
  const Int_t jptbin1[]={1,2,3,4,5, 2,3, 5,  10,11,12,15, 15};

  //old good 
  /*
  const Int_t jptbin0[]={firstjpt, 2,3,4,5};   //{firstjpt, 2,3,4,5,5,   3};//{2,3};//
  const Int_t jptbin1[]={nJpt,     2,3,4,nJpt};//{nJpt,     2,3,4,6,nJpt,8};//{2,8};//
  */
  const Int_t nbj = hn.Contains("IDFF") ? sizeof(jptbin0)/sizeof(Int_t) : 1;

  for(Int_t ieta = 1; ieta<=nEta+1; ieta++){
    TString tmptag;

    if(ieta<=nEta){
      hh->GetAxis(axEta)->SetRange(ieta, ieta);
      tmptag=Form("%s_ETA%d", tag.Data(), ieta-1);
    }
    else{
      hh->GetAxis(axEta)->SetRange(1,4);
      tmptag=Form("%s_ETAall", tag.Data());
    }

    for(Int_t ibj = 0; ibj<nbj; ibj++){
      const Int_t b0 = jptbin0[ibj];
      const Int_t b1 = jptbin1[ibj];

      jax->SetRange(b0, b1);

      OutputHist(hh, tmptag+Form("_JPT%.0f-%.0f", jax->GetBinLowEdge(jax->GetFirst()), jax->GetBinUpEdge(jax->GetLast())));
    }
  }

  printf("\n");

}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=6){
    printf("bad argc! %d\n", argc); 
    return 1;
  }

  const Int_t tpccutmode = atoi(argv[3]);
  const Int_t tofmode = atoi(argv[4]);
  const Int_t kMC  = atoi(argv[5]);

  const TString tag(Form("%s_TPC%dTOF%d_kMC%d",argv[2], tpccutmode, tofmode, kMC));

  /* http://root.cern.ch/root/html/TH1.html
     When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
     h->SetDirectory(0);          for the current histogram h   
     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference  
     When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are auto \
     matically deleted.    */
  //very important
  TH1::AddDirectory(kFALSE);

  TFile *_file0 = new TFile(argv[1]);
  printf("\nFile content:\n");
  _file0->ls();

  TString dn;
  TString ln;

  if(tofmode==0){//no TPC, TOF at all
    axPreID = -999;
  }
  else if(tofmode>=1 && tofmode<=4){//both TPC and TOF
    axPreID = axComb;
  }
  else if(tofmode==9){//only TPC
    axPreID = axtpcid;
  }
  else{
    printf("wrong tofmode %d\n", tofmode); return 1;
  }
  printf("checkcheck tofmode %d axPreID %d\n", tofmode, axPreID);

  //anyway if tofmode==0 or ==4 the TOF doesnot have an effect
  const Int_t tofid = ( (tofmode==0||tofmode==9)? 3 : tofmode);
  if(tag.Contains("kMC1")){//mc
    dn=Form("PWGJE_FragmentationFunction_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
    ln=Form("fracfunc_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
  }
  else{//data
    //gsi data
    //dn=Form("PWGJE_FragmentationFunction_jetsAOD_FASTJET04_B0_Filter00768_Cut00150_jetsAODMC2_FASTJET04_B0_Filter00768_Cut00150_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
    //ln=Form("fracfunc_jetsAOD_FASTJET04_B0_Filter00768_Cut00150_jetsAODMC2_FASTJET04_B0_Filter00768_Cut00150_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);

    //lego data
    dn=Form("PWGJE_FragmentationFunction_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_noGenJets_trackTypeUndef_jetTypeUndef_tpc%d_tof%d_cl0", tpccutmode, tofid);
    ln=Form("fracfunc_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_noGenJets_trackTypeUndef_jetTypeUndef_tpc%d_tof%d_cl0",tpccutmode, tofid);

  }
  printf("please check\n%s\n%s\n", dn.Data(), ln.Data());

  TDirectory *fdir=_file0->GetDirectory(dn);
  if(!fdir){
    printf("no fdir!!\n");exit(1);
  }
  printf("\nDirectory content:\n");
  fdir->ls();

  TList * ll=(TList*)fdir->Get(ln);
  if(!ll){
    printf("no list!!\n");exit(1);
  }
  printf("\nList content:\n");
  ll->ls();

  //THnSparseD * THnIncl=0x0;
  THnSparseD * THnIDFF=0x0;

  printf("\nGetHist ...\n");
  GetHist(ll, THnIDFF, "THnIDFF");
  //GetHist(ll, THnIncl, "THnIncl");

  THnIDFF->GetAxis(axZ)->SetRangeUser(0,1-1e-3);

  printf("\ngetOutputHist ...\n");
  getOutputHist(THnIDFF, tag);
  //getOutputHist(THnIncl, tag);

  printf("done!!\n");
  return 0;
}
